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ABSTRACT 



Context. The improvement in observational facilities requires refining the modelling of the geometrical structures of astrophysical 

objects. Nevertheless, for complex problems such as line overlap in molecules showing hyperfine structure, a detailed analysis still 

requires a large amount of computing time and thus, misinterpretation cannot be dismissed due to an undersampling of the whole 

space of parameters. 

Aims. We extend the discussion of the implementation of the Gauss-Seidel algorithm in spherical geometry and include the case of 

hyperfine line overlap. 

Methods. We first review the basics of the short characteristics method that is used to solve the radiative transfer equations. Details 

are given on the determination of the Lambda operator in spherical geometry. The Gauss-Seidel algorithm is then described and, by 

analogy to the plan-parallel case, we see how to introduce it in spherical geometry. Doing so requires some approximations in order 

to keep the algorithm competitive. Finally, line overlap effects are included. 

Results. The convergence speed of the algorithm is compared to the usual Jacobi iterative schemes. The gain in the number of 

iterations is typically factors of 2 and 4 for the two implementations made of the Gauss-Seidel algorithm. This is obtained despite the 

introduction of approximations in the algorithm. A comparison of results obtained with and without line overlaps for N2H + , HCN, 

and HNC shows that the 7 = 3-2 line intensities are significantly underestimated in models where line overlap is neglected. 



1. Introduction 

Our understanding of the evolution of the molecular complexity 
in astrophysical media relies both on the quality of the observa- 
tions and on their interpretation. The increasing spatial resolu- 
tion that will be provided by the new generation of ground- and 
space-based telescopes (Herschel, ALMA, JWST) will enable 
^^ . finer geometrical structures to be infered in the objects under 
Vh ' study. Hence, the best way to interpret the origin of the molecu- 
C3 lar line emission and to derive physical and chemical conditions 
is to use non local radiative transfer calculations with realistic 
geometries adapted to the angular resolution. While the Sobolev 
approximation can be safely used when a large velocity gradi- 
ent is present, interstellar clouds often show linewidths that are 
thermal in nature or that involve low-velocity fields, so the use 
of such an approximation is rather doubtful, and a widely used 
method in the field of molecular astrophysics is the Monte-Carlo 
technique (Bernes, 1979). This method is particularly attractive 
due to its conceptual simplicity and because it is easily imple- 
mented to treat multidimensional geometries. Nevertheless, an 
inherent drawback of the method is that it suffers from numer- 
ical noise and presents a slow convergence rate when optically 
thick lines ar e modelled. Although line overlaps can be easily 
implemented dGonzalez-Alfonso & Cernicharol ll993), comput- 
ing time becomes very expensive when treating a case such as 
N2H + where the two nitrogen atoms contribute to the hyperfine 
structure. 
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Since the 90s, various techniques have been introduced to 
solve the coupled set of radiative transfer and statistical equi- 
librium equations. These methods lead to substantial gains in 
term of calculation speed and accuracy of the solutions. Among 
them, the preconditioning of statistical equilibrium equations 
entails drastic improvement in the convergence properties of 
the iterati ve scheme for problems w ith large optical depths 
(see iRvbicki & Hummed Il99ll [l992). A crucial issue is that 
the Lambda iteration scheme does not alwa ys converge to 
the true so lution for large optical depths (see Mih alasl 119781 : 
Auer, 1991). Another progress consisted in the introduction by 
Truiillo Bueno & Fabiani Bendicho ( 1995) of the Gauss-Seidel 
algorithm for the solution of radiative transfer problems. 

The description of this method was first made for the 
case of one-dimensio nal plane-pa r allel s te llar atm ospheres 
(Truiillo Bueno & Fabiani Bendicho, 1995; Paletou & Leger, 
2007) and was further generalised to the case of two and 
three-dimensional radiative tran sfer problems with mu l tilevel 
atoms in C a rtesian coordinates (Fabiani Bendicho et al., 1997; 
iLegeretaL, [2 007 ) an d to the case of scattering line po- 
larization (Truiillo Bueno & Manso Sainz, 199 9]). I t was then 
described (Asensio Ramos & Trujillo Bueno, 2006) and used 
dCernicharo et al.. 2006) in spherical geometry. Although the ge- 
ometry of interstellar clouds cannot always be assumed to be as 
simple as spherical, it enables a first guess more appropriate than 
plan-parallel geometry in order to model centrally condensed 
clouds. In this article, we give more details on the implementa- 
tion in spherical geometry, and we extend the discussion to the 
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case of line overlap. We find that in spherical geometry, the re- 
quirements of the original Gauss-Seidel algorithm cannot be ful- 
filled without increasing the computational time (even in the case 
of isolated lines). To address this problem, we introduced several 
approximations into the algorithm. Interestingly, despite these 
approximations, the usual increase in reaching convergence in 
terms of number of iterations persists. Moreover, these approxi- 
mations do not entail false convergence, in comparison to results 
obtained within the A iteration scheme. 

Several molecules used to trace the physical conditions in 
molecular clouds harbour a hyperfine structure that c omplicates 
the a nalysis of the observed emission ( see e.g. iDaniel et all 
2007). Indeed, the influence of the overlap is usually not con- 
sidered when interpreting observations despite the fact that for 
a molecule with hyperfine structure, the line blending increases 
with the rotational quantum number. For high rotational tran- 
sitions, all the hyperfine lines are usually overlapped in astro- 
physical objects. Previous works have shown that, for HCN, line 
overlaps play an important role in the pumping of molecular 
levels and have to be considered in order to interpret the rela- 
tive intensities of hyperfine transitions ( cf. Guilloteau & Baudrv, 
I198U iGonzalez-Alfonso & Cernicharol 119931) . One goal of this 
paper is to extend the discussion to other molecules, i.e. HNC 
and N2H + , and to derive the general radiative transfer effects 
that are common to the different species. 

This article is organised as follows. In section [27X1 we first 
overview the short characteristics (SC) method. In section 12.21 
the main equations used in the preconditioning of statistical 
equations are given. In section 12.31 the implementation of the 
Gauss-Seidel algorithm in spherical geometry is discussed. In 
section [3] a test case for H2CO gives the main properties of the 
various algorithms. Finally, in section [4] the influence of intro- 
ducing line overlap in the calculations is discussed for the N2H + , 
HCN, and HNC molecules. 

2. Theory 

2. 1 . Short characteristics in spherical geometry 



Iv(V,) 





Fig. 2. Geometric quantities in spherical geometry. 



Assuming a parabolic trend for the source functions with respect 
to the opacity of the line, we can derive the intensity analytically 
as follows: 



Mi = on S v (t,-_i ) + Pi S v (t,-) + ji S v (t,+i ) 



(2) 



where the expressions of the coefficients a,, jS,-, and y, are given 
by Olson & Kunasz ( 1987) and depend on At, and At,-+i. These 
latter quantities are derived assuming that the absorption coeffi- 
cients k v (s) vary linearly between two consecutive grid points. 

In spherical geometry, the radiation field is fully charac- 
terised introducing two variables, r and 0. Numerically, one has 
to adopt a discretized space for these two variables. At a given ra- 
dius r, and for an angle #,, the specific intensity is known within 
the SC scheme from the geometrical quantities represented in 
Fig. [2]. In the case #, < n/2, we find 



d\ t i - -r,-cos0 + Jrf +l - rf sin 2 8j 
02 i - arcos 



d 2 +r 2 - r 2 ^ 
"1,1 ^ 'i+i ' i 



2d u r i+l t 
di,t = n cos Gi - yjr 2 k - r] sin 2 t 



(3) 



03j - arcos 



- 2 - rl - d 2 



with 



k = 



2 d 2 j r k 
( i+l if t = n/2 
i if 0,e]7r/2; £,,_![ 

i - 1 if 0t > fti-i 



(4) 
(5) 



Fig. 1. Transfer along a short characteristic. 



The angle £ m „ is defined according to 



• equation , we use the SC method intro- 
duced by lOlson & Kunasa(ll987l) . Here, we briefly overview the 
method and give the main equations in order to introduce quan- 
tities that will be used in the following discussion. 

If one considers the radiation propagation along a ray that 
is discretized (cf. Figure [TJ, the specific intensity b etween two 
consecutive points varies as dOlson & Kunaszj [l987l) 



I v (Ti) = / v (Ti-i)e- AT '+ f's v (T)e T <- T dT 



o /,■ = /,.ie 



■Mi 



(1) 



£ m> „ = n - arccos \\\ - \ — 

1 'm 



(6) 



and stands for the angle sustained by the layer of radius r„ seen 
from the layer of radius r m > r n . 

The quantities d\j and dij are used to derive At,- and At,-+i. 
The angle 02j is used to interpolate the incoming intensity. 
Moreover, the angle 0t,j is used in the determination of At,+i 
when a radial velocity field is introduced. For the case t > n/2, 
the above formulae also apply, using the angle n - 0j and then 
inverting respectively 02j with #3,, and d\j with d%,h 
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2.2. Accelerated Lambda Iteration 



Following Ols on et alj d!986l) . rather than the lambda iteration 
scheme (LI), we use an alternative scheme known as the accel- 
erated (or approximated) lambda iteration (ALI) since this one 
presents better convergence properties. In that scheme, the spe- 
cific intensity obeys 



Ivo = (A v n-K n )[St]+Kn\fv]+Xn 



(7) 



where A v q stands for the exact Lambda operator and A* n for an 
approximated Lambda operator (ALO). In this expression, the 
indexes correspond to the frequency v and to the direction of 
propagation Q. = (r, 0). The quantities indexed by + are deter- 
mined using the level populations obtained at the previous iter- 
ation. The value ,TJ" n corresponds to the external radiation field 
which is transmitted to the point under consideration. From Eq. 
|7]the LI scheme would be obtained letting A* n = 0. 

Several studies aiming at determining the most effective 
ALO to be used so that the convergence is reached in the fastest 
way have been conducted. Non diagonal Lambda operators were 
found to be efficient in reducing the number of iterations but 
at a computational cost per iteration higher than for a diagonal 
Lambda operator. Finally, the use of a diagonal ALO seems to 
lead to the fastest conve rgence rates (see e.g. iPuls & Herrerol 
U988tlQlsoneTaflll986h . 

When we use the diagonal part of the full lambda opera- 
tor for A* n , we find that instead of solving the usual statisti- 
cal eq uilibriu m equations ( SEE), one has to solve the system 
dRvbicki & Hummed 1 1 99 lh 

-eff. 



= £[»/Aji(l - A*,.) - imBij - njBjdj'jt ] 

j>i 

- J][n*Ay(l - Ay) - (njBji - mB^Jy] 

j<< 

+Yj [n J c J i ~ niCi j ] 



j*i 



with 



jf = J + .. -A* xSt: 



(8) 



(9) 



In the above equation, Sy is the usual frequency-independent 

source function for an isolated molecular line. Th term 7, corre- 
sponds to the radiation field averaged over frequency and an- 
gle. The above expressions apply with or without any source 
of background continuum radiation. When there is a source of 
back ground continuum radiat ion, the source function is given 
by dRvbicki & Hummed. [l99lh 



S v = r y 5 y + (l -nj)S c 
with 

Kij(v) 

r u — 

K U (V) + K C (V) 



(10) 



(ID 



where /cy(v) is the absorption coefficient due to the molecular 
transition and k c (v) is the continuum absorption coefficient. The 
averaged Lambda operator is given by 



A y = iXI 0y(v)A: ^ dvdQ ' 

where 0,y(v) is the intrinsic line shape function. 



(12) 



2.2.1. Overlap of hyperfine lines 

The overlap between hype rfine lines is treated as described by 
Rybicki & Hummer (1992). Here we give a brief overview of 
the main equations. In the case of overlapping lines, the source 
function without background continuum radiation is given by 



Sy = ——}Kij(v)Si 



(13) 



y 
>>j 



where Sy is the usual frequency independent source function 
encountered in the case of an isolated line and k v is the total 
absorption coefficient at frequency v given by 



K(V) = 2_j «ij( v ) 



(14) 



y 



Within the SC scheme, these quantities are used to perform 
the analytical integration along a characteristic using Eqs. Q] 
and [2] As explained in Rybic ki & Hummer] dl992l) . the precon- 
ditioning of the statistical equilibrium equations can be made 
in different ways. In the present work, we use the precondi- 
tioning with the same transition. This corresponds to eq. 2.23 
of Rvbicki & Hummer ( 1992). Rewriting this latter equation in 
terms of the Lambda operator rather than the T operator which 
acts on emissivities and considering that a diagonal ALO acts 
like a simple multiplication leads to the same form of the pre- 
conditioned statistical equilibrium equations as given by Eq. [8] 
The only difference is that the averaged diagonal ALO is now 
given by 



A; 



4tt J n J v 



<P U (v) 



Kjj(v) 
K(V) 



A* n dvdQ 



(15) 



We note that the full preconditioning strategy might be more 
adapted to the case of the hyperfine structure where lines are 
close in frequency. Thus, the present choice is motivated by the 
SEE equations remaining unchanged and are still given by Eq. 

i 

In the appendix, the way the Lambda operator is constructed 
is presented and the differences with the plan-parallel case are 
emphasised. 

2.3. Gauss-Seidel 

The Gauss-Seidel (GS) iterative scheme was introduced in ra - 
diative transfer by Trujillo Bueno & Fabiani Bendichol {1995). 
The superiority of this algorithm was then shown, by compar- 
ison to Jacobi-based ones, in order to rapidly reach convergence 
even for the complex case of multilevel atoms in two and three- 
dimensional Cartesian geometry (see iFabiani Bendicho et all 
1997). Recently, a detailed description of the implementation 
of the GS algorith m for multilevel atomic (or molecular) case 
was presented by iPaletou & Legerf d2007l) . So far, the GS al- 
gorithm has been extensively described in plan-parallel geom- 
etry. The descripti on of the implementation in spher ical geome- 
try has been done lAsensio Ramos & Truiillo Bue no (2006), and 
here we extend the discussion made in this article. Before, we 
will review the main attributes required by the algorithm. 

At each iteration, the grid is first swept from the outermost 
layer until the innermost one (say, from k - N — > 1 ), and 
during this process, part of the averaged radiation field is 
computed, i.e. for 6 e [0;7r/2]. Then, once the innermost 
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layer is reached, the grid is swept in the other direction. This 
leads to the computation of the averaged radiation field that 
corresponds to 6 e \tij2\ti\. Within the Jacobi iterative scheme, 
all the grid is swept until reaching the outermost layer, and 
then statistical equilibrium equations are solved. Within the GS 
iterative scheme, once the averaged radiation field is known for 
a layer, the statistical equilibrium equations are solved for it. 
At this stage, various updates have to be made to account for 
the newly derived populations. Using parabolic interpolation, 
this implies correcting both the outward and inward radiation 
field respectively, for the current and the next grid point. These 
corrections are made so that when deriving the populations for 
a layer k, the populations considered are those obtained during 
the current iteration for the layers i < k and those corresponding 
to the previous iteration for the layers i > k, noted jn" CM 'J and 

We now describe how the GS algorithm should be imple- 
mented in spherical geometry and how we have implemented it 
in practice. In what follows, we note If, the intensity 4(6*) for 
6 e [0; jt/2] and If" the intensity I k (0) for e [n/2; n]. 



2.3.1. Upward inversion 

Let us consider first that the grid has been swept from the outer- 
most layer to the innermost one. This leads to the computation of 
II f | . Once the innermost layer is reached, we can compute 

the intensities If" since the incoming intensities are given by If 
(If there is a central source of continuum radiation, the intensi- 
ties are an input parameter.). We then know J, and we can invert 
the statistical equilibrium equations for this layer. Before going 
to the point k — 2, there are three corrections that in principle 
have to be done to conform with the requirements of the original 
GS algorithm: 



which is computed with n" ew and «2 W 



and n° ld and n° ld 



- an update of I" 

- an update of I° ut 

- an update of 7™ computed with n 

Note that the first correction should be performed because the 
intensities If are used to derive the incoming intensities when 
determining 1°'" . However, this correction is not necessary in 
practice since it influences the convergence rate just slightly (cf. 
section[3]l. Once these corrections are done, we go forward in the 
propagation. 

Suppose we reach the point k where new populations are ob- 
tained. Let us now consider the corrections that are needed to 
have an exact GS iterative scheme. Since new populations have 
been obtained for the layer k, a requirement of the algorithm is 
that, at the point k+ 1, the populations should be determined with 
{n new } for j e[l;k] and {n old } for j e[k+l;N]. Now, consider 
the way If" depends on the various populations and so, through 
the dependance on the incoming intensities. For this purpose, 
we refer to Fig. [3] which shows the dependencies of intensities. 
We see that to obtain the intensities If* 1 , it is necessary to know 
the inward intensities If, as well as the outward intensities I?"' v 
which in turn depend on 7?"L etc. Finally, we see that fully ac- 
counting for the new population n'f w in the computation of If" 
would require first to perform the propagation along the inward 
direction, leading to the updated set {/"}jg[i^t+i]. Then, the prop- 
agation would be done again along the outward direction so that 
the set {I° u '}je[Uk] would effectively be given considering the new 
populations nf w . In practical applications, such a process would 



lead to computational times longer than those obtained within 
the Jacobi iterative scheme. Moreover, it should be noted that 
this problem does not arise in plan-parallel geometry: in that 
case, once the population is updated at point k, the computation 
can be performed at point k — 1 and will take exactly the new 
populations into account (cf. Fig. [3]), since the radiation field at 
the innermost boundary is an input parameter of the model. 

To make the algorithm computationally competitive, we pro- 
ceed as follows. By analogy with the plan-parallel case, a way to 
implement the Gauss-Seidel algorithm would be not to account 
for the influence of the new populations n" k ew on the incoming in- 
tensities If_ j and If when we compute the outgoing intensities 
If" . This would lead to performing 

- an update of Iff y using the new populations at points k - 
2, At— 1, and k using Iff 2 and If 1 for the incoming intensities, 
calculated with the old populations. 

- an update of If" using the populations n"ff v n 
As for the previous correction, the incoming intensities 7? 
are given by the old populations. 

- an update of If +l using the populations n 



new „nH n 0,d 
k • andB *+l- 



..new 
'k ' 



old 
'k+V 



old 
jfc+2 - 



As previously said, the first correction is in fact not necessary 
for speeding up the convergence rate. Actually, the mixing of old 
and new populations for the point k reduces slightly the conver- 
gence rate (cf. section[3]l. Finally, we only perform the two latest 
corrections, which is equivalent to the corrections described by 
Asensio Ramos & Trujillo Bueno (2006). 



2.3.2. Upward and downward inversion 

As explained above, within the usual GS scheme and for 
each iteration, the grid is first swept from k — N to k - 
1 and then from k - 1 to t = N, and during this 
second step, s tatistical equilibrium equations are solved . As 
mentioned by Trujillo Bueno & Fabiani B endichol (1 19951) and 
Trujil lo Bueno & Manso Sainzl (11999), a way to increase the 
convergence speed is to invert the statistical equilibrium equa- 
tions during both the upward and downward sweeping of the 
grid. 

During the downward propagation and when the iterative 
process has been initiated, we already know the intensities 
{If"}je[i;N]- Thus, when reaching point k, we directly compute 
the averaged radiation field for 8 e [0; n]. First we compute If 
using If 



and then 1°'" using I'," and 7°"' as incoming intensities. Here we 



'£-1 

note that the intensities If" x in fact involve a mixing of the pop- 
ulations obtained during the two last iterations. This is because 
the first correction described in the previous section is not per- 
formed and that the intensities are saved after the second correc- 
tion. Finally, since the intensities If" are computed at this stage, 
there is only one correction needed so that, after inversion, we 
perform an update of If using the populations n™f v nf/ w , and 
nf^ v During the downward propagation, we have to save the 
intensities If calculated with the new populations since these 
intensities will be used during the upward propagation. The up- 
ward sweeping of the grid is similar to the one described in the 
previous section. 



3. Numerical results 

The numerical implementation was checked against 
the results obtained with t he c ode described by 
Gonzalez-Alfonso & Cernicharo (1993). Various models 
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Spherical Geometry 
I k (0); 0£[V2;tt] 

— > I k (0) ; 0e[O;7r/2] 
L^ I k ^(0); 0£[tt/2;7t] 

— > Ij^l(fl); 0£[O;tt/2] 
— > I k _ 2 (0); 0£[tt/2;tt] 

— > I k _ a (0); 0e[O;7T/2] 



I 2 (0); 0£[O;tt/2] 



Plan paralell Geometry 
I k (0); 0e[TT/2;TT] 

■♦ I k ^(0); 0£[tt/2;tt] 

— > I k _ a (fl); 0£[tt/2;tt] 

— > I k _ 3 (0); G [7t/2;tt] 



IjCe); 0e[rr/2;^] 

I > 1^0); 0£[O;tt/2] 



li(0); 0e[7T/2;Tr] = I ext 



Fig. 3. Dependence of the intensities. 



were used in the comparison, covering different physical 
conditions, as well as low, moderate, and large optical depths. 
The differences in the derived level populations ranged from 
a tenth of a percent to a few percent when integrated optical 
depths were higher than ~ 100. This corresponds to the expected 
accuracy of numerical rad iative transfer codes as explained in 
Ivan Zadelhoff et al.1 (120021) . 

To accelerate the convergence we used the ALI method de- 
scribed in Sect. 12.21 . In spherical geometry, the diagonal of the 
full Lambda operator is given by Eqs. IB. 101 IB. Ill and IB. 121 
However, determining the exact diagonal is time consuming: 
first, it requires the evaluation of the opacities At* + , that de- 
pend on all the internal shells crossed by the rays and second, 
it requires interpolating in angle the coefficients c* , to obtain 



them at the values 6* , 



n+l 



In practice, we thus chose not to keep 
the corresponding terms and we approximated the diagonal of 
the Lambda Operator just keeping, for 9 > n/2, 



■At; 



- if m > K + 1 : 

Ki.m - /? m (0m)+7 m -l(0~-i)e" 

- if m — K : 

A m , m = a m (0 m ) +P m (6Q + \Pm(.0 + m ) + JmifQ] e 



-ATt,, 



(16) 



(17) 



For 9 < it 12, we conserve the terms given by expressions IB. II 
and EI] 

These terms permit us to accelerate the rate of convergence 
notably by comparison to the LI scheme. Moreover, the approx- 
imation made in deriving the diagonal of the ALO does not in- 
duce false convergence. Both aspects can be seen in Fig. |4] where 
the maximum difference with respect to the converged solution 
is plotted. The model considered here corresponds to a sphere 
of uniform physical parameters with n(H;j) = 10 4 cm" 3 and tem- 
perature T = 25K. The sampling consists in ~ 500 spatial grid 
points and 20 angles. This has been found to insure a conver- 
gence for the level populations better than 0. 1 % for the first 12 
levels. The molecule considered is ortho-F^CO, and 14 energy 



levels were included in the calculation. The molecular parame- 
ters were taken from the references given in Table Q] Since no 
analytical solution exists for this problem, we used as reference 
populations those obtained with the LI scheme. The error rela- 
tive to the converged solution for the considered g rid, noted C e , 
is defined according to Eq. 20 of lAuer et al.l d!994l) . 

Figure |4] also shows the rate of convergence of the different 
algorithms in terms of the iteration number. For the GS algo- 
rithms, two curves are shown: the grey one corresponds to the 
case where the incoming intensities are corrected (i.e. correction 
1 of Sect. |2.3.Tb and the black one corresponds to the case where 
this correction is omitted. We see that, for this model, the ALI 
scheme speeds up the convergence by a factor 1 .3 in comparison 
to the LI scheme. The GS down and GS up-down enabled us to 
reduce the number of iterations by factors of 2.0 and 4.4 respec- 
tively by comparison to the ALI scheme. In Sect. 12. 31 the various 
corrections to be done for these two algorithms were discussed. 
For each correction, the code calls a formal solver that takes the 
source functions of the lines as input and outputs the averaged 
radiation field and Lambda operator, as well as the intensities. 
Thus, considering the number of calls to the formal solver and 
in the different schemes, we have for each layer 2 calls in the ALI 
scheme, 4 calls in the GS down scheme, and 6 calls in the GS 
up-down scheme. The number of calls gives the relative time per 
iteration of each algorithm. In practice, the GS schemes imple- 
mentation can be improved by performing some bookkeeping, 
in particular for the incoming intensities. These lead to times per 
iteration, in the GS down and GS up-down schemes, which are 
respectively 1.45 and 2.3 times greater than in the ALI scheme. 
These times were corrected for the time spent in the code in in- 
ternal storage on the hard drive. Such a procedure is done in our 
implementation of the GS schemes while not for the ALI one 
and makes the GS scheme ~ 20% slower than permitted by the 
algorithmic. Moreover, these longer times are compensated for 
by the reduced number of iterations required (see above). Thus, 
for the current implementation, the GS down and GS up-down 
respectively reduce the total calculation time by factors of 1.4 
and 1.9. The numerical implementation of the GS scheme could 
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molecule 


rate coefficients 


dipole moment 


energy levels 


H,CO 
HCN 
HNC 
N 2 H + 


Green (1991) 

Monteiro & Stutzki (1986) 

Monteiro & Stutzki (1986) 

Daniel et al. (2005) 


Fabricant et al. (1977) 

Ebenstein & Muenter (1984) 

Blackman et al. (1976) 

Havenith et al. (1990) 


Miiller et al. (2005) 
Miiller et al. (2005) 
Miiller et al. (2005) 
Casellietal. (1995) 



Table 1. References for the molecular parameters used in this work. 



still be improved to gain some computing time. The corrections 
in spherical geometry are formall y similar to the one in p lan- 
parallel geometry, and for this case. iPaletou & Legerl d2007l) find 
that the time per iteration of the GS down scheme was only 1.3 
longer than that of the ALI scheme. Our code has been written 
in a very modular and structured way. By decreasing the num- 
ber of internal calls to routines and functions, we could still gain 
computing time and approach the 1.3 value of plan-parallel ge- 
ometry. 

Nevertheless, the GS schemes become more advantageous 
if, rather than parabolic interpolation, we use linear interpola- 
tion. In that case, there are fewer corrections to carry out so that 
the time per iteration of the GS down and GS up-down scheme 
are only 1 .4 and 1 .7 that of the ALI scheme. The drawback is that 
if using linear interpolation, one needs to use more grid points 
in order to converge to the true solution. Nevertheless, for the 
model considered here, we find that the differences in level pop- 
ulations are inferior to 1%. In Fig. [5] the line parameters of the 
most opaque lines are plotted for the two types of interpolations 
and the results are indistinguishable on the scale of the graph. 
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Fig. 5. Excitation temperature versus opacity for a model con- 
sidering ortho-H2CO. 
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Fig. 4. Convergence properties of the various algorithms in terms 
of number of iterations. 



We note that the various algorithms discussed here can still 
be improv ed if combined with the Ng acceleration technique 
(Ng, 1974). For the various models tested, we found that it typ- 
ically leads to a reduction by a factor of 3 in the number of it- 
erations. Al ternatively, the successive over-relaxation technique 
can be used Trujillo Bueno & Fabiani Bendicho ( 1995). 



4. Overlap of hyperfine lines 

The numerical implementation of the overlap 

was checked according to the results reported by 

iGonzalez- Alfonso & Cernicharo! (Il993l) . For the molecules 



considered in what follows the molecular parameters are sum- 
marised in Table [TJ Figure [6] shows the variation in excitation 
temperature versus line opacity for the three HCN hyperfine 
lines associated with the 7=1-0 rotational transition. 
The physical parameters of the mo dels were taken from 
IGonzalez- Alfonso & Cernicharo! ( 1 19931) . The line parameters 
are given with (solid lines) and without (dashed lines) intro- 
duction of the over lap. By comparison to the resu l ts rep orted 
in Figs. la-Id of IGonzalez- Alfonso & Cernicharo! d 1993b . we 
see that the qualitative behaviour of T ex is correctly reproduced 
for all the hyperfine lines when the overlap is introduced. 
Nevertheless, a quantitative comparison of the line parameters, 
with and without overlap, shows differences in the order of 5— 
10%. These differences are due to the choice of grid parameters, 
especially the number of radial grid points, adopted to reach 
convergence. 

Figure|7]shows the influence of the overlap on the modeliza- 
tion of the hyperfine lines associated with the 3 lowest N2FL 
rotational transitions. This molecule is particularly interesting 
for infering the physical conditions of the densest parts of cold 
dark clouds. Indeed, it was found that this molecule does not suf- 
fer from strong depletion for Fb volume densities below ~ 10 6 
cm" 3 . The model p arameters are identical to those reported in 



Dani el et all d2007) for the source L63. In this figure, we see 
that the 7=1-0 line remains mainly unchanged with differ- 
ences in the line fluxes lower than 5%. The differences increase 
when considering higher rotational lines. For the 7 = 3-2 line, 
including the overlap increases the line flux by ~ 20-30 %. This 
can be understood as follows. For a group of blended hyperfine 
lines, the average radiation field seen by each line is greater when 
considering the overlap, since more molecules can interact with 
the radiation. Consequently, the pumping of the upper level of 
each transition is more efficient, which results in an enhancement 
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Fig. 6 . Line parameters obtained for t he models described in 
iGonzalez- Alfonso & Cernicharo! d 1993b 
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Fig. 7. Effect of the overlap on the line profile of L63. The black 
line correspond to a model without overlap and the red lines cor- 
respond to models that include overlap. 



of the populations of the highest rotational levels. Therefore the 
overla p changes the physical parameters reported in Danie letall 
(2007). Mainly, H2 volume densities and/or temperatures were 
overestimated in this study. 

The effect induced by line overlap obviously depends on in- 
dividual hyperfine line opacities, as well as on the frequency off- 
set between the lines. Figure [8] shows the emerging intensities 
for the three lowest rotational lines of HCN, HNC, and N 2 H + . 
The models consist of a uniform density core of diameter 1'. 
The density is n(H2) = 10 5 cm 3 , the temperature T=10 K, and 
the turbulence is set to v rur t, = 0.15 km.s -1 . For each molecule, 
calculations are performed at abundances of 10~ 10 , 2.5 ■ 10~ 9 , 
and 2.5 • 10~ 8 . First, it can be seen that, at the lowest abundance, 
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Fig. 9. Effect of overlap on relative line intensities. The black 
curve corresponds to a model that includes line overlap and the 
red line to a model without overlap, the latter being scaled by the 
integrated intensity ratio. 



the effect introduced by the overlap is weak for N2I-F and HCN. 
The effect is stronger for HNC due to closer hyperfine frequen- 
cies. Second, we find that the effect introduced by the overlap 
is approximatively proportional to the column density, as long 
as individual line opacities < 1 . In the present calculations, this 
is found for the 7 = 3-2 lines of each species for the 2 low- 
est abundances. For these lines, we see that the integrated in- 
tensity ratios of the calculations with and without overlap are 
approximatively proportional to the abundance ratio. In a more 
optically thick regime, the increase in the integrated intensity 
ratio depends less on the column density and finally reaches a 
saturation regime where the intensity ratio starts to be affected 
by self-absorption effects. This can be seen in Fig. [8] for the 
7 = 2-1 and 7=1-0 lines. Finally, the error introduced in the 
interpretation of observations when neglecting the line overlap 
can be especially important for HNC. Indeed, we see in Fig. [8] 
that, for conditions typical of dark clouds, including line over- 
lap leads to an intensity increase of up to a factor 2. Moreover, 
for similar abundances, the overlap differentially affects HCN 
and HNC and thus, neglecting it may lead to overestimating the 
X(HNC)/X(HCN) ratio. 

Another effect introduced by the overlap is to modify the 
relative intensities of the hyperfine components of a given ro- 
tational transition. Figure [9] shows the 7 = 2-1 line of N2H + 
for a model with abundance 2.5 10~ 9 . The line obtained with- 
out overlap is represented multiplied by a factor of 1.54, which 
corresponds to the integrated intensity ratio obtained previously. 
In this figure the position in frequency of the different hyperfine 
components is also represented as are their line strengths. These 
are given in percent by reference to the value obtained by sum- 
mation over all the hyperfine components. Considering the red 
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carried out did not show false convergence by comparison 
with LI calculations. 

The implementation of the Gauss-Seidel algorithm in spher- 
ical geometry also requires some approximations. These are 
introduced by analogy with the plan-parallel case. 
The convergence rate of the GS down and GS up-down are 
similar to what is obtained in plan-parallel geometry, i.e. 
they reduce the number of iterations to reach convergence 
by factors of 2 and 4 respectively. The gain in terms of com- 
puting time strongly depends on the numerical implementa- 
tion, and in the present code, the GS up-down reduces the 
calculation time by a factor ~ 2 by comparison to the ALI 
scheme. 

Line overlaps can lead to a large increase in the predicted 
intensities, the effect being greater for higher rotational lines. 
Especially, differential effects for HCN and HNC can alter 
the estimate of the abundance ratio X(HNC)/X(HCN), and 
models where line overlap is neglected will overestimate it. 



satellite, we see that, when including the overlap, the intensity 
associated with the strongest transition (group 1 at V ~ 2.5 km 
s -1 ) is reduced by comparison to the group of weaker lines sit- 
uated at v ~ 3 km s _1 (group 2). This effect on the profile for 
the N2H 4 " 7 — 2—1 red sate llite agrees with the observations 
reported in Daniel et al. (2007). Two effects yield to the current 
emerging profile. The main effect is introduced by the number of 
strong transitions that overlap in each group. In group 1, the line 
strengths of the two weakest transitions are low enough so that 
there is mostly one line that contributes to the emission. On the 
other hand, group 2 corresponds to 4 lines with equivalent line 
strengths. Thus, for this group of lines, the pumping of the upper 
levels is enhanced producing an increase in the emergent inten- 
sities. Another effect of the overlap is to create a flow of popula- 
tion that depends on the relative line strengths of the overlapping 
lines. Figure[10]shows excitation temperature ratios between the 
transitions of group 1 obtained with (solid lines) and without 
(dashed lines) introduction of the overlap. We see in this figure 
that the overlap enhances the excitation temperatures of the 2 
weakest lines. This means that the levels JF\F = 223 and JF\F 
= 221 are more efficiently populated when considering the over- 
lap due to the presence of the JF\F - 212- 101 transition that 
is close in frequency. This can be understood when considering 
the way the average radiation field is obtained when solving the 
radiative transfer. For the strongest line, the presence of weak 
hyperfine components close in frequency does not have a strong 
influence since the opacities/intensities in the wings of the in- 
trinsic line profiles are low. The average radiation field is thus 
only slightly modified. On the other hand, for a weak line that is 
close to the frequency of a strong line, the average radiation field 
is increased. This will pump the corresponding energy levels and 
result in an increase in intensity for the weakest lines. 

5. Conclusions 

We have presented technical details concerning the determina- 
tion of the Lambda operator. Moreover, we gave details on the 
implementation of the Gauss-Seidel algorithm in spherical ge- 
ometry. Finally, the numerical code was applied to the treatment 
of line overlap. The main conclusions are: 

1 . The diagonal of the Lambda operator, as obtained within the 
short characteristics method, has to be approximated in prac- 
tical applications. The use of approximated diagonal terms 
still enable speeding the convergence. Moreover, the tests 
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Appendix A: Discretized Lambda operator: 
plan-parallel geometry 

As pointed out by Rybicki & Hummer] dl991l) the Lambda op- 
erator is in fact an approximation of the true Lambda operator 
since it depends on the discretization of the structure of the ob- 
ject. Thus, the Lambda operator is a square matrix whose dimen- 
sion corresponds to the number of spatial grid points. 

Within the SC scheme, the Lambda operator can be con- 
structed easily. As a first step, from Eqs.[T]and|2]one finds that 
along a characteristic 



h = / e T <+ Y ; A/,-*e 



fc=0 



(A.l) 



A' 



where Ar,j = r,- - T/. To obtain a compact expression, we adopt 
from now on the notation aj l = »,-, cq = /?,-, and a\ - y,- for the 
terms that appear in Eq. [2] Introducing it in the previous expres- 
sion, we find 



i+i ( 1 



^n^^ 



with 



k=o V=-i 



Ki = 



Sv(T k ) 



if 0<k + l<i 



u k+l 

otherwise 



(A.2) 



(A.3) 



In the above expression, the source functions are given at optical 
depths that correspond to the intersection of the ray with the dis- 
cretized structure of the cloud. In order to construct the Lambda 
operator, it is then necessary to introduce a notation where the 
source functions are given as a function of grid points. 



a=n 



Vffl 



i,-,(e) . 



Fig. A.l. Notation adopted in plan-parallel geometry. 

In plan-parallel geometry, the values for S v (t>) are readily 
assigned to values of the source functions for each layer. From 
Fig. IA.ll we have in the case 9 e [0; n/2] that the point k — in 
expression |A.2| corresponds to the layer n — 1 . The point k = i— 1 
corresponds to the layer n - m. In the above expression, we thus 
first introduce the variable n = k + 1 and then a notation where 
all the quantities are indexed according to the layer, i.e. 
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U n-\+l 
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(A.4) 



The conditionlA.3ltransforms as 



^ = 



2 n+l 



if 1 <n + I <m+ 1 



otherwise 



(A.5) 



where the coefficient a , keeps the same definition given above, 
but is now labelled according to the layer index. It leads to 

m+l I 1 



n=\ 



V/=-l 



s: 



(A.6) 



A„ 



In this expression, A m „ stands for the element of the m' h line and 
n' h column of the Lambda operator when e [0; n/2]. 

Equivalent transformations can be done in the case 8 e 
[tt/2; 7r]. Defining N so that n = N — k, identifying i — N ' — m and 
introducing the transformation 
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leads to 



N ( 1 



a H! =Z Z^ e 



S n v . (A.9) 

n=m—l V/=-l > 

This final expression defines the element A,„,„ when 6 e [n/2; n]. 



Appendix B: Discretized Lambda Operator: 
spherical geometry 

In spherical geometry, the same procedure of the plan-parallel 
case can be followed to obtain the components of the Lambda 
operator. In this case, however, all the shells are not necessarily 
crossed for a given direction, and furthermore, some shells are 
crossed two times before reaching the current point. 




Fig. B.l. Definition of geometrical quantities. 

We consider a cloud discretized in Af concentric shells. In 
the case 9 e [0;7r/2], the transformations are similar to the sec- 
ond case examined in the previous section with two exceptions. 
First, in the SC frame, the coefficients b +l , depend on an angle 
defined in a coordinate system associated to the layer n + l. This 
dependence was not outlined in the plan-parallel case since all 



F. Daniel and J. Cernicharo: Solving radiative transfer with line overlaps using Gauss Seidel algorithms. 1 1 

the coefficients along a characteristic are at a unique angle value, - if k — K : 

namely 9. In what follows, we introduce the notation 0* +J , which _ . + r . _ 

refers to the angle to be taken into account in the evaluation of A »>,k = &m,K (#) + Pk{9 + k)z t " ,k + [(^k(&k) + Pk(@k)) e T " 

b + [ , (see Fig. IB. II ). The superscript in this notation distinguishes . //i+\ a -AT + „l ., /i ^ Q /■ w 

acute (-) and obtuse (+) angles. J 

Secondly, depending on the value of 9, the source function x^ quantity j[ mk (6 /) [ s defined according to 
corresponding to the point i + 1 can be associated either to shell 
m or to shell m— 1. In this case, we obtain for 9 e [0; n - £ m , m _i] 
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In the above expressions, we have introduced the notation 



In the expressions derived above, the dependence on the vari- 



Ar m,«- A r ay that intersects the shells m and n defines two seg- a bl e 9 is not made explicit but would appear in the definitions of 
ments. The opacity along the longest one defines Ar+ „, and the q± +[ an d A T f n k+r The expressions obtained for A„ a are suitable 
opacity along the smallest one defines Ar ffl „ (see Figs. EI] and for construct ing an operator of whatever bandwidth. However, 

we are interested here in the diagonal of this operator: 




Fig. B.2. Definition of geometrical quantities. 



Now, let us consider the case 9 > n/2, and designate K the 
innermost shell that is intersected by the ray. This shell is such 
that we have i; m x < 9 < g m ,K-i, with the angle £„,# defined 
according to Eq. [6] 

Manipulations of Eq. |A~2l lead to 
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ABSTRACT 



Context. To investigate the physical nature of the 'nucleated instability' of proto giant planets, the stability of layers in static, radiative 

gas spheres is analysed on the basis of Baker's standard one-zone model. 

Aims. It is shown that stability depends only upon the equations of state, the opacities and the local thermodynamic state in the 

layer. Stability and instability can therefore be expressed in the form of stability equations of state which are universal for a given 

composition. 

Methods. The stability equations of state are calculated for solar composition and are displayed in the domain -14 < lgp/[gcnr 3 ] < 

0, 8.8 < lge/[ergg~'] < 17.7. These displays may be used to determine the one-zone stability of layers in stellar or planetary 

structure models by directly reading off the value of the stability equations for the thermodynamic state of these layers, specified by 

state quantities as density p, temperature T or specific internal energy e. Regions of instability in the (p, e)-plane are described and 

related to the underlying microphysical processes. 

Results. Vibrational instability is found to be a common phenomenon at temperatures lower than the second He ionisation zone. The 

/(--mechanism is widespread under 'cool' conditions. 



Key words, giant planet formation - ^--mechanism - stability of gas sphe 



1. Introduction 

In the nucleated instability (also called core instability) hypothe- 
sis of giant planet formation, a critical mass for static core enve- 
lope protoplanets has been found. Mizuno ( 1980) determined the 
critical mass of the core to be about 12 M e (M e = 5.975 X 10 27 g 
is the Earth mass), which is independent of the outer boundary 
conditions and therefore independent of the location in the solar 
nebula. This critical value for the core mass corresponds closely 
to the cores of today's giant planets. 

Although no hydrodynamical study has been available many 
workers conjectured that a collapse or rapid contraction will en- 
sue after accumulating the critical mass. The main motivation 
for this article is to investigate the stability of the static envelope 
at the critical mass. With this aim the local, linear stability of 
static radiative gas spheres is investigated on the basis of Baker's 
( 1966) standard one-zone model. 

Phenomena similar to the ones described above for giant 
planet formation have been found in hydrodynamical mod- 
els concerning star formation where protostellar cores ex- 
plode (Tscharnuter 1987 Balluch 1988), whereas earlier studies 
found quasi-steady collapse flows. The similarities in the (mi- 
cro)physics, i.e., constitutive relations of protostellar cores and 
protogiant planets serve as a further motivation for this study. 



2. Baker's standard one-zone model 



In this section the one-zone model of Baker (1966), originally 
used to study the Cephei'd pulsation mechanism, will be briefly 
reviewed. The resulting stability criteria will be rewritten in 
terms of local state variables, local timescales and constitutive 
relations. 

Baker ( 1 19661 ) investigates the stability of thin layers in self- 
gravitating, spherical gas clouds with the following properties: 

- hydrostatic equilibrium, 

- thermal equilibrium, 

- energy transport by grey radiation diffusion. 

For the one-zone-model Baker obtains necessary conditions for 
dynamical, secular and vibrational (orpulsational) stability (Eqs. 
(34a, b,c) in Baker 1966). Using Baker's notation: 



M, 
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mass internal to the radius r 

m mass of the zone 

ro unperturbed zone radius 

po unperturbed density in the zone 

Tq unperturbed temperature in the zone 

LrQ unperturbed luminosity 

£ t h thermal energy of the zone 

and with the definitions of the local cooling time (see Fig.[T]i 

L r0 



(1) 
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Fig. 1. Adiabatic exponent Fi. H is plotted as a function of lg internal energy [ergg -1 ] and lg density [gem -3 ]. 



and the local free-fall time 



Table 1. Opacity sources. 



^ff 



3tt 4t!t 3 



/ 32G 3M r ' 

Baker's K and cr have the following form: 

n 1 
o-q = — — 
V8 T ff 

K = ^32 1 x ff 

n 6 t co ' 

where E& ~ m(Po/po) has been used and 
' — (*&), 



(2) 


Source 


77 [K] 




Yorke 1979, Yorke 1980a 

Kriigel 1971 

Cox & Stewart 1969 


< 1700" 

1700 < T < 5000 

5000 < 


(3) 

(4) 


a This is footnote a 





(5) 



is a thermodynamical quantity which is of order 1 and equal to 1 
for nonreacting mixtures of classical perfect gases. The phys- 
ical meaning of <jq and K is clearly visible in the equations 
above. <x represents a frequency of the order one per free-fall 
time. K is proportional to the ratio of the free-fall time and the 
cooling time. Substituting into Baker's criteria, using thermody- 
namic identities and definitions of thermodynamic quantities, 



V ari = 



<91nP\ 
<91np/ s 

d\nT 
dlnP 



Xp 



Xt = 



d\nP\ 
dlnp) 7 

dlnP 
d\nT 



dlriK 
dlnP 



dlnK 



d\nT )j 



one obtains, after some pages of algebra, the conditions for sta- 
bility given below: 



^1^-4) > ° 

8 T Z ff 



TcoT; 



-riV ad 



CO • ff 



1 - 3/4* p 



X T 



(k t - 4) + K P + 1 



-r?v ad 



4V ad - (Vad/Cj. + K p ) - — - 

31 1 



> 

> 



(6) 

(7) 
(8) 



For a physical discussion of the stability criteria see Baker 
( TT966T > or Cox ( fT980l 

We observe that these criteria for dynamical, secular and vi- 
brational stability, respectively, can be factorized into 

1. a factor containing local timescales only, 

2. a factor containing only constitutive relations and their 
derivatives. 

The first factors, depending on only timescales, are positive by 
definition. The signs of the left hand sides of the inequalities ([6j, 
© and (O therefore depend exclusively on the second factors 
containing the constitutive relations. Since they depend only on 
state variables, the stability criteria themselves are functions of 
the thermodynamic state in the local zone. The one-zone stabil- 
ity can therefore be determined from a simple equation of state, 
given for example, as a function of density and temperature. 
Once the microphysics, i.e. the thermodynamics and opacities 
(see Table[TJ), are specified (in practice by specifying a chemical 



composition) the one-zone stability can be inferred if the ther- 
modynamic state is specified. The zone - or in other words the 
layer - will be stable or unstable in whatever object it is imbed- 
ded as long as it satisfies the one-zone-model assumptions. Only 
the specific growth rates (depending upon the time scales) will 
be different for layers in different objects. 

We will now write down the sign (and therefore stability) 
determining parts of the left-hand sides of the inequalities ©, 
© and dSj an d thereby obtain stability equations of state. 

The sign determining part of inequality © is 3Fi - 4 and it 
reduces to the criterion for dynamical stability 



Ft>-. 

Stability of the thermodynamical equilibrium demands 

X p >0, c„ >0, 

and 

Xt >0 

holds for a wide range of physical situations. With 

r 3 - 1 = l*L > o 

pT c v 



v a 



i 



ad 



r, 



> o 



(9) 

(10) 

(11) 

(12) 
(13) 
(14) 



we find the sign determining terms in inequalities (0 and ([8]) 
respectively and obtain the following form of the criteria for dy- 
namical, secular and vibrational stability, respectively: 



3r, 



* dyn 



> o 



1 - 3/4^ 
X T 



(k t - 4) + K P + 1 =: S sec > 



4V ad - (V ad /c r + k p ) - — =: S vib > . 
jl i 



(15) 
(16) 

(17) 



The constitutive relations are to be evaluated for the unperturbed 
thermodynamic state (say (p ( >, Tq)) of the zone. We see that the 
one-zone stability of the layer depends only on the constitu- 
tive relations Ti, V ad , % T , x p , « p , k t . These depend only on the 
unperturbed thermodynamical state of the layer. Therefore the 
above relations define the one-zone-stability equations of state 
S dyn , Ssec and S v ib- See Fig.|2]for a picture of 5 v ib- Regions of 
secular instability are listed in Table 1 . 
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Fig. 2. Vibrational stability equation of state S v ib(lg e, lgp). > 
means vibrational stability. 

3. Conclusions 

1 . The conditions for the stability of static, radiative layers in 
gas spheres, as described by Baker's (1966) standard one- 
zone model, can be expressed as stability equations of state. 
These stability equations of state depend only on the local 
thermodynamic state of the layer. 

2. If the constitutive relations - equations of state and 
Rosseland mean opacities - are specified, the stability equa- 
tions of state can be evaluated without specifying properties 
of the layer. 

3. For solar composition gas the k- mechanism is working in the 
regions of the ice and dust features in the opacities, the H2 
dissociation and the combined H, first He ionization zone, 
as indicated by vibrational instability. These regions of in- 
stability are much larger in extent and degree of instability 
than the second He ionization zone that drives the Cephei'd 
pulsations. 
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